﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace Nature
{

    public abstract class ModelBoundPotentialConservationLaw : Nature.Numerics.ConservationLowBase
    {
        protected readonly ModelObject Model;
        public readonly int Nsp;

        protected ModelBoundPotentialConservationLaw(ModelObject model)
        {
            this.Model = model;
            this.Nsp = model.Nsp;
        }        
    }

    public abstract class IdealGasPotentialConservationLaw : ModelBoundPotentialConservationLaw
    {
        protected abstract double CalcValue(double variable, double[] composition);

        public readonly IIdealGasSpeciesThermoPotentials Thermo;

        private double[] m_composition;

        protected double[] SpeciesPotentials { get; private set; }

        public virtual void Configure(double fixedValue, ReadOnlyArray<double> composition)
        {
            base.FixedValue = fixedValue;
            composition.CopyTo(m_composition);
        }

        public void Configure(double fixedValue, ReadOnlyArray<decimal> composition)
        {
            base.FixedValue = fixedValue;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_composition[isp] = (double)composition[isp];
            }
        }

        protected IdealGasPotentialConservationLaw(ModelObject model)
            : base(model)
        {
            this.Thermo = model.IdealGasSpeciesThermoPotentials;
            m_composition = new double[Nsp];
            SpeciesPotentials = new double[Nsp];
        }

        protected override double CalcValue(double variable)
        {
            return this.CalcValue(variable, m_composition);
        }
    }

    public abstract class IdealGasPotentialConservationLawSolver : Nature.Numerics.Basic.Newton<IdealGasPotentialConservationLaw>
    {
        protected IdealGasPotentialConservationLawSolver(IdealGasPotentialConservationLaw conservationLaw)
            : base(conservationLaw) { }

        public static IdealGasPotentialConservationLawSolver CreateSpecificInternalEnergyConservationLawTemperatureSolver(ModelObject model)
        {
            return new IdealGasSpecificInternalEnergyConservationLawTemperatureSolver(model);
        }

        public static IdealGasPotentialConservationLawSolver CreateMolarInternalEnergyConservationLawTemperatureSolver(ModelObject model)
        {
            return new IdealGasMolarInternalEnergyConservationLawTemperatureSolver(model);
        }

        public static IdealGasPotentialConservationLawSolver CreateSpecificEnthalpyConservationLawTemperatureSolver(ModelObject model)
        {
            return new IdealGasSpecificEnthalpyConservationLawTemperatureSolver(model);
        }

        public static IdealGasPotentialConservationLawSolver CreateMolarEnthalpyConservationLawTemperatureSolver(ModelObject model)
        {
            return new IdealGasMolarEnthalpyConservationLawTemperatureSolver(model);
        }

        public static Nature.Numerics.Basic.Newton<IdealGasMolarEntropyConservationLawTemperatureEquation> CreateMolarEntropyConservationLawTemperatureSolver(ModelObject model)
        {
            return new IdealGasMolarEntropyConservationLawTemperatureSolver(model);
        }
    }

    class IdealGasSpecificInternalEnergyConservationLawTemperatureEquation : IdealGasPotentialConservationLaw
    {
        public IdealGasSpecificInternalEnergyConservationLawTemperatureEquation(ModelObject model) : base(model)
        {}

        protected override double CalcValue(double variable, double[] yk)
        {
            double usum = 0d;
            var uk = SpeciesPotentials;
            Thermo.SpecificU(variable, uk);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                usum += uk[isp] + yk[isp];
            }
            return usum;
        }
    }

    class IdealGasMolarInternalEnergyConservationLawTemperatureEquation : IdealGasPotentialConservationLaw
    {
        public IdealGasMolarInternalEnergyConservationLawTemperatureEquation(ModelObject model)
            : base(model)
        { }

        protected override double CalcValue(double variable, double[] xk)
        {
            double usum = 0d;
            var uk = SpeciesPotentials;
            Thermo.MolarU(variable, uk);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                usum += uk[isp] + xk[isp];
            }
            return usum;
        }
    }

    class IdealGasSpecificEnthalpyConservationLawTemperatureEquation : IdealGasPotentialConservationLaw
    {
        public IdealGasSpecificEnthalpyConservationLawTemperatureEquation(ModelObject model)
            : base(model)
        { }

        protected override double CalcValue(double variable, double[] yk)
        {
            double hsum = 0d;
            var hk = SpeciesPotentials;
            Thermo.SpecificH(variable, hk);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                hsum += hk[isp] + yk[isp];
            }
            return hsum;
        }
    }

    class IdealGasMolarEnthalpyConservationLawTemperatureEquation : IdealGasPotentialConservationLaw
    {
        public IdealGasMolarEnthalpyConservationLawTemperatureEquation(ModelObject model)
            : base(model)
        { }

        protected override double CalcValue(double variable, double[] xk)
        {
            double hsum = 0d;
            var hk = SpeciesPotentials;
            Thermo.MolarH(variable, hk);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                hsum += hk[isp] + xk[isp];
            }
            return hsum;
        }
    }

    public class IdealGasMolarEntropyConservationLawTemperatureEquation : IdealGasPotentialConservationLaw
    {
        double[] m_logxk;

        public IdealGasMolarEntropyConservationLawTemperatureEquation(ModelObject model)
            : base(model)
        {
            Pressure = Constants.Atmosphere;
            m_logxk = new double[Nsp];
        }

        public override void Configure(double molarEntropy, ReadOnlyArray<double> composition)
        {
            Configure(molarEntropy, Constants.Atmosphere, composition);
        }

        public void Configure(double molarEntropy, double pressure, ReadOnlyArray<double> composition)
        {
            base.Configure(molarEntropy, composition);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_logxk[isp] = Math.Log(1.0e-30 + composition[isp]);
            }
            Pressure = pressure;
            LogPressureAtm = Math.Log(Pressure / Constants.Atmosphere);
        }

        public double Pressure
        {
            get;
            private set;
        }

        public double LogPressureAtm
        {
            get;
            private set;
        }

        protected override double CalcValue(double variable, double[] xk)
        {
            double ssum = 0d;
            var sk = SpeciesPotentials;
            Thermo.S(variable, sk);
            for (int isp = 0; isp < Nsp; ++isp)
            {
                ssum += xk[isp] * (sk[isp] - m_logxk[isp]);
            }
            return (ssum - LogPressureAtm) * Constants.Rgas;
        }        
    }

    public class IdealGasIsentropicCompressionTemperatureEquation : Nature.Numerics.IFunction
    {
        readonly ModelObject Model;
        readonly IIdealGasSpeciesThermoPotentials Thermo;
        readonly int Nsp;

        double[] m_logxk;
        double[] m_sk;
        double[] m_xk;
        double m_compressionRatio;
        double m_temperature0;
        double m_ssum0;

        public IdealGasIsentropicCompressionTemperatureEquation(ModelObject model)
        {
            Model = model;
            Thermo = model.IdealGasSpeciesThermoPotentials;
            Nsp = model.Nsp;
            m_logxk = new double[Nsp];
            m_sk = new double[Nsp];
            m_xk = new double[Nsp];
            m_compressionRatio = 1.0d;
            m_temperature0 = 298.15;
        }

        public void Configure(double molarEntropy, double temperature0, double compressionRatio, ReadOnlyArray<double> x)
        {
            m_compressionRatio = compressionRatio;
            m_ssum0 = molarEntropy / Constants.Rgas;
            m_temperature0 = temperature0;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_xk[isp] = x[isp];
                m_logxk[isp] = Math.Log(1.0e-30 + m_xk[isp]);
            }
        }

        #region IFunction Members

        public double Compute(double temperature)
        {
            Thermo.S(temperature, m_sk);
            double ssum = 0d;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                ssum += m_xk[isp] * (m_sk[isp] - m_logxk[isp]);
            }
            double func = ssum - m_ssum0 - Math.Log(m_compressionRatio * temperature / m_temperature0);
            return func;
        }

        #endregion
    }

    public class IdealGasMolarEntropyConservationLawPressureCalculator
    {
        public readonly ModelObject Model;
        public readonly IIdealGasSpeciesThermoPotentials Thermo;
        public readonly int Nsp;

        double[] m_logxk;
        double[] m_sk;

        public IdealGasMolarEntropyConservationLawPressureCalculator(ModelObject model)
        {
            Model = model;
            Thermo = model.IdealGasSpeciesThermoPotentials;
            Nsp = model.Nsp;
            m_logxk = new double[Nsp];
            m_sk = new double[Nsp];            
        }

        public double CalcPressure(double temperature, double pressure, double molarEntropy, ReadOnlyArray<double> x)
        {
            for (int isp = 0; isp < Nsp; ++isp)
            {
                m_logxk[isp] = Math.Log(1.0e-30 + x[isp]);
            }
            return CalcPressure(temperature, pressure, molarEntropy, x, m_logxk);
        }

        public double CalcPressure(double temperature, double pressure, double molarEntropy, ReadOnlyArray<double> x, ReadOnlyArray<double> logx)
        {
            Thermo.S(temperature, m_sk);
            double ssum = 0d;
            for (int isp = 0; isp < Nsp; ++isp)
            {
                ssum += x[isp] * (m_sk[isp] - logx[isp]);
            }
            double ssum0 = molarEntropy / Constants.Rgas + Math.Log(pressure / Constants.Atmosphere);
            return pressure * Math.Exp(ssum - ssum0);
        }
        
    }

    class IdealGasSpecificInternalEnergyConservationLawTemperatureSolver : IdealGasPotentialConservationLawSolver
    {
        public IdealGasSpecificInternalEnergyConservationLawTemperatureSolver(ModelObject model) 
            : base(new IdealGasSpecificInternalEnergyConservationLawTemperatureEquation(model))
        { }
    }

    class IdealGasMolarInternalEnergyConservationLawTemperatureSolver : IdealGasPotentialConservationLawSolver
    {
        public IdealGasMolarInternalEnergyConservationLawTemperatureSolver(ModelObject model)
            : base(new IdealGasMolarInternalEnergyConservationLawTemperatureEquation(model))
        { }
    }

    class IdealGasSpecificEnthalpyConservationLawTemperatureSolver : IdealGasPotentialConservationLawSolver
    {
        public IdealGasSpecificEnthalpyConservationLawTemperatureSolver(ModelObject model)
            : base(new IdealGasSpecificEnthalpyConservationLawTemperatureEquation(model))
        { }
    }

    class IdealGasMolarEnthalpyConservationLawTemperatureSolver : IdealGasPotentialConservationLawSolver
    {
        public IdealGasMolarEnthalpyConservationLawTemperatureSolver(ModelObject model)
            : base(new IdealGasMolarEnthalpyConservationLawTemperatureEquation(model))
        { }
    }

    class IdealGasMolarEntropyConservationLawTemperatureSolver : Nature.Numerics.Basic.Newton<IdealGasMolarEntropyConservationLawTemperatureEquation>
    {
        public IdealGasMolarEntropyConservationLawTemperatureSolver(ModelObject model)
            : base(new IdealGasMolarEntropyConservationLawTemperatureEquation(model))
        { }
    }

    public class IdealGasIsentropicCompressionTemperatureSolver : Nature.Numerics.Basic.Newton<IdealGasIsentropicCompressionTemperatureEquation>
    {
        public IdealGasIsentropicCompressionTemperatureSolver(ModelObject model)
            : base(new IdealGasIsentropicCompressionTemperatureEquation(model))
        { }
    }
}
